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1) Introduction 

We study the optimal control of a differential linear system 

(1) = Ay + Bv, 

where the state variable y(t) belongs to H n and the control function v (•) takes its values 
in H m , with n and m beeing given integers. Matrix A is composed by n lines and n 
columns and matrix B contains n lines and m columns. Both matrices A and B are 
independent of time. With the ordinary differential equation © is associated an initial 
condition 

(2) y(0) = y 

with yo given in IR™ and the solution of system © © is parametrized by the function v (•): 
The control problem consists of finding the minimum «(•) of some quadratic functional 
J(.): 

(3) J(«(.)) < J( V (.)), V«(.). 

The functional J(») depends on the control variable function is defined by the 

horizon T > 0, the symmetric semi-definite positive n by n constant matrix Q and the 
symmetric definite positive m by m constant matrix R. We set classically : 

(4) J(v(.)) = \ £{Qy(t),y{t)) dt +l £{Rv{t),v{t)) dt ■ 

• Problem (PP) ([2D d3J) d3J) is a classical linear quadratic mathematical modelling of dy- 
namical systems in automatics (see e.g. Lewis (Le86j). When the control function v(m) 
is supposed to be square integrable (v(») G L 2 (]0; T[, H m )) then the control problem 
©©©([ID has a unique solution u(») G L 2 (]0; T[, H m ) (see for instance Lions jLl68] l 
When there is no constraint on the control variable the minimum «(•) of the functional 
J(v) is characterized by the condition: 

(5) dJ(u).w = 0, Vw G L 2 (]0, T[, H m ) , 
which is not obvious to compute directly. 

• When we introduce the differential equation © as a constraint between y(m) and v (•), 
the associated Lagrange multiplyer p(») is a function of time and is classically named the 
adjoint variable. Research of a minimum for J(.) (condition ©) can be rewritten in the 
form of research of a saddle point and the evolution equation for the adjoint variable is 
classical (see e.g. Lewis (Le86j): 

(6) ^ + A t p + Qy = 0, 
with a final condition at t — T, 

(7) p(T) = 

and the optimal control in terms of the adjoint state p(») takes the form: 

(8) Ru(t) + B t p(t) = 0. 
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• We observe that the differential system ([I])© together with the initial condition (T2|) 
and the final condition (J2J) is coupled through the optimality condition (jHJ). In practice, we 
need a linear feedback function of the state variable y(t) instead of the adjoint variable 
p(t). Because adjoint state p(») depends linearily on state variable y(») we can set: 

p{t) = X{T-t).y{t), 0<t<T, 

with a symmetric n by n matrix X(m) which is positive definite. The final condition ([7]) 
is realized for each value y(T), then we have the following condition: 

(9) X(0) = 0. 

We set K = BR^ 1 B^; we remark that matrix K is symmetric positive definite, we replace 
the control u(t) by its value obtained in relation ([8]) and we deduce after elementary 
algebra the evolution equation for the transition matrix X(m): 

dX 

(10) — - (XA + A X) + XKX - Q = 0, 

which defines the Riccati equation associated with the control problem ([TJ (j2J) (J3J) (|3j) . 

• In this paper we study the numerical approximation of differential system (I9"j) (ll(jp . 
Recall that the given matrices Q and K are n x n symmetric matrices, with Q semi- 
definite positive and K positive definite; the matrix A is an n by n matrix without any 
other condition and the unknown matrix X(t) is symmetric. We have the following 
property (see e.g. Lewis |Le86] ). 

Proposition 1. The solution of Riccati equation is positive definite. 

Let K, Q, A be given n x n matrices with K, Q symmetric, Q positive and K definite 
positive. Let X(») be the solution of the Riccati differential equation ffTU]) with initial 
condition (Q. Then X(t) is well defined for any t > 0, is symmetric and for each t > 0, 
X(t) is definite positive and tends to a definite positive matrix as t tends to infinity: 
X(t) — > Xqo if t — > oo. Matrix X^, is the unique positive symmetric matrix which is 
solution of the so-called algebraic Riccati equation: 

— (XA + A t X) + XKX - Q = 0. 

• As a consequence of this proposition it is usefull to simplify the feedback command 
law ([8]) by the associated limit command obtained by taking t — > oo, that is: 

(11) v(t) = —R" 1 B^ X OQ y(t) , 

and the differential system ([T|) (fTTj) is stable (see e.g. |Le86]). The practical computation 
of matrix I M by direct methods is not obvious and we refer e.g. to Laub [ La79] . If 
we wish to compute directly a numerical solution of instationnary Riccati equation (fTUj) 
classical methods for ordinary differential equations like e.g. the forward Euler method 

^(X J+1 - XA + Xj K Xj - (A t X j + XjA) - Q = 0, 
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or Runge Kutta method fail to maintain positivity of the iterate Xj+i at the order 
U ■ 1): 

(12) (X j+1 x,x)>0, Vi G E", x^O, 

if Xj is positive definite and if time step At > is not small enough (see e.g. Dieci and 
Eirola [DE96J). Morever, there is to our best knowledge no simple way to determine a 
priori if time step At > is compatible or not with condition (lT2~j) . 

• In the following, we propose a method for numerical integration of Riccati equation 
(TTO]) which maintains condition (fl2"|) for each time step At > 0. We present in second 
section the simple case of scalar Riccati equation and present the numerical scheme and 
its principal properties of the general case in section 3. We describe several numerical 
experiments in section 4. 



2) Scalar Riccati equation 

• When the unknown is a scalar variable, we write Riccati equation in the following 
form: 

(13) — + kx 2 - 2ax - q = 0, 

Lit- 

with 

(14) k > 0, q > 0, 
and an initial condition: 

(15) x(0) = d, d > 0. 

We approach the ordinary differential equation (13) with a finite difference scheme of 
the type proposed by Baraille |Ba91J for hypersonic chemical kinetics and independently 
with the "family method" proposed by Cariolle |Ca79] and studied by Miellou [Mi84j . We 
suppose that time step At is strictly positive. The idea is to write the approximation 
Xj + \ at time step (j + l)At as a rational fraction of xj with positive coefficients. We 
decompose first the real number a into positive and negative parts : a = a + — a~ ; 
a + = max(0; a) > 0, a~ = max(0; —a) > 0, a + a~ = and factorize the product x 2 into 
the very simple form: 

\ x ) j+i/2 = x i x i+ l ■ 
Definition 1. Numerical scheme in the scalar case. 

For resolution of the scalar differential equation (13), we define our numerical scheme by 
the following relation: 

(16) J+ At — ~ ^ ^ X 3 X 3 +1 ~ 2a + Xj + 2a~Xj + i — q = 0. 

• The scheme (TIB"]) is implicit because some linear equation has to be solved to compute 
Xj + \ when Xj is supposed to be given. In the case of our scheme this equation is linear 
and the solution Xj + i is obtained from scheme (I16p by the homographic relation: 

( 1 + 2a + At ) xj + q At 



(17) x 



j'+i 



kAtXj + (1 + 2a~ At) 
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Proposition 2. Algebraic properties of the scalar homographic scheme. 

Let [xj)j^js be the sequence defined by initial condition : Xq = x(0) = d and recurrence 
relation ( I17p . Then sequence {xj)j^ is globally defined and remains positive for each 
time step: Xj > 0, Vj G IN, VAt > 0. If At > is chosen such that: 

(18) 1 + 2|a|At - kqAt 2 ^ 0, 

then (xj)j e ]N converges towards the positive solution x* of the "algebraic Riccati equation" 

kx 2 — lax — q = 

and 

(19) x* = -(a + + hq) ■ 

• In the exceptional case where At > is chosen such that ( TLB"]) is not satisfied, then 

the sequence (xj)je]N is equal to the constant 1+ l A + t At for j > 1 and the scheme f[T7|) 
cannot be used for the approximation of Riccati equation (TT3j) . 

Theorem 1. Convergence of the scalar scheme. 

We suppose that the data k, a, q of Riccati equation satisfy (THj) and f lT8|) and that the 
datum d of condition f JT5|) is relatively closed to x* , i.e.: 



(20) - — + 77 < rf-x* < C, 

where C is some given strictly positive constant (C > 0) , x* calculated according to 
relation ( |T9|) is the limit in time of the Riccati equation, r is defined from data k, a, q 
by: i 

and 77 is some constant chosen such that 

(21) " <ri< h- 

• We denote by x(t; d) the solution of differential equation (TT5|) with initial condition 

f fT5|) . Let (xj(At; c?A))(jGiN) be the solution of the numerical scheme defined at the 
relation f TT7|) and let d& be the initial condition: 

x (At ; d A ) = g?a . 

We suppose that the numerical initial condition g?a > satisfies a condition analogous 

to flU]): 1 

-■ r- + ^ < rf A - x* < C, 
k r 

with C and i] > equal to the constant introduced in (120]) and satisfying (12T|) . 

• Then the approximated value (xj(At ; c?A))jeiN is arbitrarily closed to the exact value 
x(jAt ; d) for each j as At — > and d& — > d . More precisely, if a ^ we have the 
following estimate for the error at time equal to jAt : 

\x(jAt; d) - Xj(At; d A ) \< A (At + | d - d A |), V j G IN , < At < B 
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with constants A > 0, B > , depending on data k, a, q, r] but independent on time step 
At > and iteration j . 

• If a = , the scheme is second order accurate in the following sense: 

\x(jAt;d) - Xj(At;d A )\ <A(At 2 + | d - d A |), V j G IN , < At < B 

with constants A et B independent on time step At and iteration j. 

A direct application of the Lax theorem for numerical scheme associated to ordinary dif- 
ferential equations is not straightforward because both Riccati equation and the numerical 
scheme are nonlinear. Our proof is detailed in [DS2k|. 



3) Matrix Riccati equation 

In order to define a numerical scheme to solve the Riccati differential equation (10) with 
initial condition ([9]) we first introduce a strictly positive real number, which is chosen 
positive in such a way that the real matrix [/il — (A + A )] is definite positive: 

(22) ~(px,x) - (Ax,x)>0, Vx ^ 0. 

Then we introduce the definite positive matrix M wich depends on fi and matrix A: 

M = -ul - A. 
2 p 

The numerical scheme is then defined by analogy with relation (16). We have the following 
decomposition : 

(23) A = A + - AT 

with A + = |^I, A~ = M, a, > 0, M positive definite. Taking as an explicit part the 
positive contribution A + of the decomposition f l2"3"|) of matrix A and in the implicit part 
the negative contribution A~ = M of the decomposition ( |23|) . we get 

(24) | Xt {Xj+1 - Xj) + \^ KX ^ + x i^ KX s) + 

\ + (M t X j+1 + X j+1 M) = fiXj + Q . 

The numerical solution given by the scheme Xj + i at time step j + 1 is then defined as 
a solution of Lyapunov matrix equation with matrix X as unknown: 

S^j X -\- X Sj = Yj 

with 

(25) S 3 = l -I+^-KX J+ AtM 
and 

(26) Yj = Xj + fiAtXj + AtQ. 

We notice that Sj is a (non necessarily symmetric) positive matrix and that Yj is a 
symmetric definite positive matrix if it is the case for Xj. 
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Definition 2. Symmetric matrices. 

Let n be an integer greater or equal to 1. We define by <S n (IR), (respectively <S+(IR), 
<S+*(1R)) the linear space (respectively the closed cone, the open cone) of symmetric- 
matrices (respectively symmetric positive and symmetric definite positive matrices). The 
following inclusions <S+*(IR) C <S+(IR) C «S n (R) are natural. 

Proposition 3. Property of the Lyapunov equation. 

Let S be a matrix which is not necessary symmetric, such that the associated quadratic 
form: ]R n 3 x i — > (x, Sx) G H, is strictly positive i.e. 

S + S l e <S+*(IR). 

Then the application (p defined by : 

(27) 5 n (M) 3 X i— > <p(X) = S t X + XS G S n QR) , 

is a one to one bijective application on the space «S n (IR) of real symmetric matrices of 
order n. Morever, if matrix <f(X) is positive (respectively definite positive) then the 
matrix X is also positive (respectively definite positive): 

if <p(X) G 5+(M) , then X G 5+(IR) . 

• The numerical scheme has been written as an equation with unknown X = Xj + % 
which takes the form: (fj{X) = Yj with (fj given by a relation of the type fl27jl with the 
help of matrix Sj defined in (125|) and a datum matrix Yj defined by relation (1261) . Then 
we have the following propositions. 

Proposition 4. Homographic scheme computes a definite positive matrix. 

The matrix Xj defined by numerical scheme (1241) with the initial condition X = is 
positive for each time step At > : 

XjeS+QR), Vj>1. 

If there exists some integer m such that X m belongs to the open cone <S+*(IR) , then 
matrix X m+ j belongs to the open cone <S^*(IR) for each j. 

Proposition 5. Monotonicity. 

Under the condition i i 

-(KX^ + X^K) < (fi + — )/, 

the scheme (124|) is monotone and we have more precisely : 

(28) (o < Xj < Xoo) (o < X 3 < X j+1 < Xoo). 
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4) First numerical experiments 
4-1 Square root function 

• The first example studied is the resolution of the equation : 

(29) ^— + X 2 - Q = 0, X(0) = 

at 

with n — 2, A — 0, K = I and matrix Q equal to 

i -i w i o w i i 
i i y V o ioo J V -i i 

• We have tested our numerical scheme for fixed value At = 1/100 and different values 
of parameter u, : /i — 0.1, 10~ 6 , 10 +6 . For small values of parameter //, the behaviour 
of the scheme does not change between /j, = 0.1 and /i = 10 -6 . Figures 1 to 4 show 
the evolution with time of the eigenvalues of matrix Xj and the convergence is achieved 
to the square root of matrix Q. For large value of parameter \x (pL = 10 +6 ), we loose 
completely consistency of the scheme (see figures 5 and 6). 

fitst eigenvalue second eigenvalue 



Q = - 

* 2 





Figures 1 and 2. Square root function test. 

Two first eigenvalues of numerical solution (/i = 0.1). 



fitst eigenvalue 



second eigenvalue 





Figures 3 and 4. Square root function test. 

Two first eigenvalues of numerical solution (/i = 10~ 6 ) 
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first eigenvalue 



second eigenvalue 




100 23 Q no 400 -j'j SO 



JOB JDQ 300 400 



Figures 5 and 6. Square root function test. 

Two first eigenvalues of numerical solution (/i = 10 +6 ). 



4-2 Harmonic oscillator 

• The second exemple is the classical harmonic oscillator. Dynamical system y(t) is 
governed by the second order differential equation with command v (t) : 



(31) 



d 2 y(t) n . dy(t) 2 . . , . , 
+ 25^^ + u 2 y(t) = bv(t) 



(32) 



Y 



dt 2 ' "" dt 

This equation is written as a first order system of differential equations : 

y(t) \ dY ( o 1 \ _ / o 

dfl J ' dF = U- 2 -2 6 ) m + U«W 

In this case, we have tested the stability of the scheme for fixed value of parameter 
fi(fi = 0.1) and different values of time step At and coefficients of matrix R inside the 
cost function of relation @: 

R= ( a ° 
V a 

• We have chosen three sets of parameters : a = At = 1/100 (reference experiment, 
figures 7 and 8), a = 10 -6 , At = 1/100 (very small value for a, figures 9 and 10) and 
a = 1/100, At = 100 (too large value for time step, figures 11 and 12). Note that 
for the last set of parameters, classical explicit schemes fail to give any answer. As in 
previous test case, we have represented the two eigenvalues of discrete matrix solution 
Xj as time is increasing. On reference experiment (figures 7 and 8), we have convergence 
of the solution to the solution of algebraic Riccati equation. If control parameter a is 
chosen too small, the first eigenvalue of Riccati matrix oscillates during the first time 
steps but reach finally the correct values of limit matrix, the solution of algebraic Riccati 
equation. If time step is too large, we still have stability but we loose also monotonicity. 
Nevertheless, convergence is achieved as in previous case. 
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fitst eigenvalue 



second eigenvalue 




120 160 




40 H 120 



Figures 7 and 8. Harmonic oscillator. 

Two first eigenvalues of numerical solution (/i = 0.1, a = 0.01, At = 0.01). 



fitst eigenvalue 



second eigenvalue 




12d 160 




12(1 160 



Figures 9 and 10. Harmonic oscillator. 

Two first eigenvalues of numerical solution (// = 0.1, a = 10~ 6 , At = 0.01). 



fitst eigenvalue 



second eigenvalue 





10 31 3D 40 



Tfl 30 90 



Figures 11 and 12. Harmonic oscillator. 

Two first eigenvalues of numerical solution (// = 0.1, a = 0.01, At = 100). 
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5) Conclusion 

We have proposed a numerical scheme for the resolution of the matrix Riccati equation. 
The scheme is implicit, unconditionnaly stable, needs to use only one scalar parameter 
and to solve a linear system of equations for each time step. This scheme is convergent 
in the scalar case and has good monotonicity properties in the matrix case. Our first 
numerical experiments show stability and robustness when various parameters have large 
variations. Situations where classical explicit schemes fail to give a solution compatible 
with the property that solution of Riccati equation is a definite positive matrix have been 
computed. We expect to prove convergence in the matrix case and we will present in 
[DS2k| experiments on realistic test models such as a string of vehicles and the discretized 
wave equation. 
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